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1  A  numerical  solution  to  the  generalized  Burgers’  radial  wave  equation  has  been  developed;  it  allows  one  to 
calculate  stepwise  the  harmonic  content  of  a  finite  amplitude  wave  in  the  frequency  domain  for  the  case  of  plane, 
cylindrical,  or  spherical  geometry.  The  finite  amplitude  wave  may  have  any  initial  harmonic  content  with  arbitrary 
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domains.  A  listing  of  the  computer  program  is  included,  u 
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A  FORTRAN  COMPUTER  PROGRAM  FOR  CALCULATING 
THE  PROPAGATION  OF  PLANE,  CYLINDRICAL,  OR  SPHERICAL 
FINITE  AMPLITUDE  WAVES 

INTRODUCTION 

Several  adequate  theories  [1,2,3],  based  on  approximations  to  the 
nonlinear  wave  equation,  have  been  developed  to  describe  the  behavior  of  a 
one-dimensional  wave  of  moderate  amplitude  as  it  propagates  through  a  nonlinear 
fluid.  Most  of  these  theories,  however,  are  not  conveniently  applied  to  the 
problem  of  describing  the  propagation  when  the  wave  is  of  arbitrary  initial 
waveform  or  when  the  linear  absorption  has  an  arbitrary  frequency  dependence. 

To  handle  these  more  general  case-',  investigators  [4,5,6]  have  adapted  the 
phenomenological  model  of  Fox  and  Wallace  [7]  to  use  a  high-speed  computer 
to  calculate  the  propagation  stepwise.  The  distance  of  propagation  in  this 
model  is  divided  into  small  intervals.  The  wave  is  first  allowed  to  distort 
over  one  interval  and  is  then  corrected  to  account  for  absorption  and  geometrical 
spreading.  The  procedure  is  then  repeated  for  the  new  waveform  over  the  next 
interval.  The  use  of  small  intervals  preserves  the  interaction  between  the 
distortion,  the  absorption,  and  the  geometrical  spreading  mechanisms.  Since 
the  distortion  mechanism  is  applied  in  the  particle  velocity  domain  with 
absorption  and  geometrical  spreading  being  applied  in  the  frequency  domain, 
it  is  necessary  to  switch  back  and  forth  between  the  two  domains  during  each 
step.  Even  with  the  use  of  the  Fast  Fourier  Transform  (FFT),  this  procedure 
is  a  time-consuming  process.  In  addition,  one  must  take  special  care  in 
applying  the  distortion  mechanism  when  the  waveform  has  a  very  steep  shock-like 
portion. 

We  describe  in  this  report  a  new  procedure  for  calculating  the 
propagation  of  plane,  cylindrical,  or  spherical  finite  amplitude  waves. 

This  procedure  performs  the  stepwise  calculations  entirely  in  the  frequency 
domain,  thus  avoiding  both  the  use  of  the  FFT  and  the  steep  waveform  problems. 

We  also  describe  a  FORTRAN  computer  program  called  FAW  that  implements  the 
new  procedure. 

The  theory  behind  the  procedure  is  described  in  Sec.  I.  A  discussion 
of  the  relative  importance  of  linear  absorption  and  nonlinear  effects  is 
presented  in  Sec.  II.  In  Sec.  Ill  the  numerical  implementation  of  the  procedure 
Manuscript  submitted  November  24,  1980. 
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is  developed.  This  is  followed  in  Sec.  IV  b>  a  description  of  the  computer 
program  FAW.  Included  are  discussions  of  the  significant  FORTRAN  variable 
names,  the  major  computational  blocks,  parameter  input,  printed  output, 
number  of  harmonics  required  in  the  calculation,  and  initial  step  size. 
Section  V  contains  a  comparison  of  results  obtained  using  the  new  procedure 
with  results  obtained  using  the  phenomenological  model.  The  report  concludes 
with  appendices  containing  sample  output  and  a  program  listing  of  FAW. 

I.  THEORY 

An  approximate  nonlinear  wave  equation  valid  for  one-dimensional  and 
progressive  plane,  cylindrical,  and  spherical  waves  in  a  lossless  fluid 
is  the  generalized  Burgers'  equation  [8] 


-I 


f  +  -  bD  -  o. 


(1) 


where  U  =  particle  velocity, 

r  =  spatial  coordinate, 

x  =  retarded  time  (t  ), 

c  ’ 
o 

cq  =  small  signal  sound  speed, 
2 


b  =  6/c 


B 


8  “  °  +  2A  )»  tJle  nonlinearity  parameter,  and 

a  =  0  (plane  waves), 

=  1/2  (cylindrical  waves),  or 
=  1  (spherical  waves). 


Equation  (1)  is  modified  to  include  linear  absorption  by  noting  that  in  the 
absence  of  nonlinearity  the  amplitude  decays  according  to 


U(r)  =  Uq  (rQ/r)  exp  [-  a(u)(r-rQ)] 


or 


=  -(a/r)  U  -  ce(w)  U. 


The  first  term  on  the  right-hand  side  of  Fq.  (3)  represents  the  loss  due  to 
geometrical  spreading  and  is  already  included  in  Fq.  (1).  The.  second  term 
represents  the  loss  due  to  linear  absorption.  Adding  this  term  to  Fq.  (1) 
results  in  the  following  nonlinear  equation  for  a  lossy  fluid  with  arbitrary 
frequency-dependent  absorption  where  dispersion  has  been  neglected: 


|y  +  (a/r)U  -  bU  (3U/3t)  =  -  ct(co)  V, 


('*) 


We  now  choose  as  a  trial  solution  a  Fourier  series  of  linear  damped  waves 
of  arbitrary  phase  with  amplitudes  that  are  a  function  of  the  spatial 
coordinate  r: 


U(r,T)  =  ^  (r^/r)3  |  sinfku^t)  +  cosfk^t)  j.  exp[-  a^fr-r^],  (5) 

1=1 


til 

where  is  the  absorption  coefficient  appropriate  for  the  k  harmonic. 

The  fundamental  frequency  u  is  chosen  less  than  or  equal  to  1/t  ,  where  t 

o  o  o 

is  either  the  period  of  the  initial  waveform  at  r=rQ  when  the  waveform  is 
periodic  or  it  is  a  time  length  sufficiently  long  to  contain  the  resulting 
waveform  at  all  desired  distances  when  the  waveform  is  transient.  Substitution 
of  Eq.  (5)  into  Eq.  (4)  yields  two  coupled  differential  equations  governing 
the  behavior  of  the  amplitude  components  G^  and  as  a  function  of  the  spatial 
coordinate  r: 
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t—  sin  (kw  T)exp  (  — « 
dr  o  K 


(r-r  )] 
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and 


O  3Hk 

2-^cos  (kv}  exp  t~°k(r'ro)1 


k=l 


Dili  r  A 

=  — -  (r  /r)a  /  .  m  exp  (-(a  »+<*  )(r-r  )] 

2  O  <.  m  o 

£,m 

|(HPG  +  GpH  )  cos  [(£-hn)co  t] 

(  -C  m  -C  m  o 

+  (HPG  -  GPH  )  cosl(£-m)w  t]i  . 

•C  in  -c  m  o( 

Factoring  out  terms  of  the  same  frequency  (i.e.,  terms  in  which  £+ro  =  k, 
£-m  =  k,  and  £-m  =  -k)  results  in 


(7) 
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k-1 

IT m(C.  C  -H.  H  )  exp  I~(a 
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+a  -a  )(r-r  )] 
k-m  m  k  o 


+  Y 'n  (C.  c  +  It  .11  )  exp  [-(a.  ■ **  >(r-r  )] 

/  t  k+m  n  k+m  m  k-hn  m  k  o 

m=l 

00 

m  (G  .  G  +  I!  .II  )  exp  [-(«  -«.)(r-r  )]( 

/  j  m-k  m  m-k  m  m-k  m  k  o  I 

m=k+l 


(P) 


and 


3H,  b(« 

k  _  o 

~3r  ,  2~ 


/k- 1 

(r„/r)aE-  <"k-A  +  Ck-»V  '-"i/.l*'".11 

lm=l 


+  >  m(IL  G  -  G.  1i  )  exp  [-(«,  -a.  )(r-r  )] 

/  -i  Vhn  m  k+m  m  *  k+m  m  k  o 


(9) 


m=l 


CO 

+  m(H  ,  G  -  G  ,11  )  exp  {-(«  .-¥*■  -«.)(r-r  )]>. 

/  j  m-k  m  ra-k  m  ’  m-k  m  k'  o  1 
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The  first  sum  on  the  right-hand  side  of  Fqs.  (8)  and  (9)  represents  contributions 

to  the  ktk  harmonic  due  to  the  interaction  of  lower  harmonics  producing  a  sum 

frequency  component  at  the  ktk  harmonic.  The  second  sum  in  the  equations 

corresponds  to  the  interaction  of  higher  harmonics  producing  a  difference 

frequency  component  at  the  ktk  harmonic.  The  last  sum  in  the  equations  repre- 

t  h 

sants  the  loss  to  the  k  harmonic  due  to  its  interaction  with  all  of  the  harmonics. 
Combining  the  last  two  sums  in  Fqs.  (P)  and  (9)  yields 


3G.  bu)  I 

1?  "  IT  (ro/r)a|2^m  (Ck-mCm  ‘  \-mV  '^k.Vk^-V’ 


m=l 


-k  >  ((W\n  +  "khnV  «■» 


on) 


k*fm  m  k' 


m=l 


and 


31!  bo) 

m=l 


m  O'  C  +  H  G,  )  exp  [-(“  +“  -«*  )(r-r  )1 

k-m  m  m  k-m  1  k-m  m  k  o 


Co 

+k  /  (H  G  -  11,  G  )  exp  [-(a,  +a  -a  )(r-r  )]  (  . 
Z—rf  m  k+m  k+n  m  k+m  m  k  o  f 

m=l  ) 


Ol) 


Equations  ( lu)  and  (11)  are  the  coupled  nonlinear  equations  that  are  numerically 
integrated  to  obtain  t 
spatial  coordinate  r. 


integrated  to  obtain  the  harmonic  amplitudes  G^  and  1'^  as  a  function  of  the 


II.  RELATIVE  IMPORTANCE  OF  LINEAR  ABSORPTION  AND  NONLINEAR  FFFFCTS 
The  Goldberg  number  r  is  defined  as 


r=ar 


(12) 
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where  a  is  the  absorption  coefficient  and  t  is  the  discontinuity  distance, 
the  point  at  which  the  waveform  would  shock  if  linear  absorption  were  absent. 
The  discontinuity  distance  is  geometry  dependent  and  for  initially  sinusoidal 
waves  is  given  by 


£p=W£T  wave), 


2  4 

c  c 

l‘c  =  ro  +  WZT  +  7-1V2  (cylindrical  wave), 
00  4r  S  U  w 
o  00 


l  =  r  exp  (c  /£U  w  r  )  (spherical  wave). 

S  O  O  OOO 


If  the  Goldberg  number  is  greater  than  unity,  the  nonlinear  effect  becomes 
important  and  shocks  are  likely.  In  this  case,  as  the  waveform  approaches  the 
discontinuity  distance,  nonlinear  effects  dominate  the  loss  due  to  linear 
absorption.  The  amplitudes  of  all  harmonics  above  the  fundamental  increase  at 
the  expense  of  the  fundamental.  After  the  discontinuity  distance  is  reached, 
however,  linear  absorption  plays  an  increasingly  larger  role  and  eventually  the 
amplitudes  of  all  the  harmonics  decrease  with  distance. 

The  use  of  the  Goldberg  number  is  important  in  deciding  on  the  number  of 
harmonics  to  retain  in  the  calculation.  If  the  Goldberg  number  is  small  compared 
to  unity,  then  the  nonlinear  effect  is  small  and  the  waveform  is  not  going  to 
shock.  In  this  case,  a  small  number  of  harmonics  will  adequately  describe  the 
waveform  at  any  position.  However,  if  the  Goldberg  number  is  large,  then  shocks 
are  likely  and  a  large  number  of  harmonics  must  be  retained  in  the  calculation. 

A  discussion  of  the  relative  error  associated  with  the  number  of  harmonics  retained 
in  the  calculation  is  found  in  Section  IV. E. 
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III.  NUMERICAL  IMPLEMENTATION 


A.  Series  Truncation 

In  order  to  numerically  integrate  Eqs.  (10)  and  (11)  the  infinite 
series  on  the  right-hand  side  of  each  equation  must  be  truncated  in  such  a 
manner  that  no  instability  is  introduced  into  the  algorithm.  This  problem  can 
best  be  examined  by  assuming  that  the  phase  of  the  initial  waveform  is  such 
that  all  the  amplitude  coefficients  1?^  are  zero  and  that  i  harmonics  are  retained 
in  the  calculation.  This  results  in  Fq.  (11)  vanishing  and  Eq.  (10)  reducing 
to 


3C.  hw 
k 


^r 


hw  0=4 

-  (r  / r)a</  m  C.  C.  exp  [-(a.  -fa  -a,)(r-r  )1 
2  o  k-m  m  1  k-m  m  k  o 

'  m=  1 

jk  i 

_k  XlCk-HnCm  6XP  f-<0k-Hn‘hlm-0k)(r-ro)]}- 

m=  1  * 


(1A) 


The  simple  truncation  used  in  obtaining  Fq.  (IA)  is  insufficient  when 
an  attempt  is  made  to  examine  the  propagation  of  the  waveform  beyond  the  dis¬ 
continuity  distance.  In  calculating  the  propagation,  the  flow  of  energy  from 
lower  to  higher  harmonics  stops  with  the  last  harmonic  retained  in  the  series. 

This  is  obvious  from  the  fact  that  the  second  series  in  Eq.  (1A)  vanishes  for 
the  harmonic.  Thus  the  use  of  simple  truncation  eliminates  the  primary 
nonlinear  energy-loss  mechanism  of  the  last  (j^')  harmonic.  The  harmonics 
preceding  the  last  are  affected  in  a  similar,  hut  less  severe,  manner.  The 
calculated  values  for  tue  last  harmonics  become  abnormally  large  relative  to 
the  lower  harmonics.  Being  "toe  large”,  these  harmonics  then  cause  an  abnormal 
growth  of  the  next  lower  harmonics  so  that  eventually  even  the  lowest  harmonics 
are  significantly  in  error.  This  instability  is  circumvented  in  the  program  by 
artificially  increasing  the  loss  of  the  last  few  fand  least  significant)  harmonics 
by  requiring  that  their  amplitude  never  exceed  the  amplitude  of  the  next  lower 
harmonic. 
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B.  Integration  Method 


Equations  (10)  and  (11)  are  stepwise  numerically  integrated  by  the 
first-order  Runge-Kutta  method  to  obtain  the  amplitudes  of  the  harmonics  at 
progressively  increasing  distances.  This  method,  as  applied  to  the  numerical 
solution  of  the  problem 

°k'  -  -Tf  -  Vr-  n-  ,f>  n7) 

and 

3H 

-~fr  -  \<r,  C,  V)  a n 


yields  the  amplitudes  at  the  K+l^1  step  as 

Ck(?:+1)  =  Ck(K)  +  [c'(M)  +  O'OM-1)]  (3°) 

and 

I'k(K+l)  -  I!r(K)  +  |  |>'(K)  +  U'(K+1)],  (?0) 

where 

O'(N)  =  Rk  fr(N),  G(N),  H(K)  \  (?1) 

!''(::)  =  SR  frfK),  0(H),  !*(K)  1,  (??) 

G'(K+1)  =  Rj^  fr(K+l),  G(N)  +  bo'(K),  »'(N)  +  hI?'(K)\  (??) 

Hk(N+l)  =  Sk  (r(f'+l),  GOO  +  hG'(K),  »'00  +  hl’"(K)*,  (?4) 
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and  h  is  the  incremental  step  size  1 9 ) .  This  procedure  has  the  advantage  of 
not  requiring  the  calculation  of  any  derivatives  of  and  as  would  be 
necessary  in  a  Taylor-series  expansion.  The  procedure  also  allows  the  step 
size  to  be  easily  changed  at  any  point  in  the  calculation.  The  only  disadvantage 
is  that  the  right-hand  sides  of  Eqs.  (17)  and  (18)  must  be  evaluated  twice 
at  each  step. 

C.  Variable  Step  Size 

In  calculating  the  propagation  of  a  finite  amplitude  wave  it  is  desirable 
to  use  the  largest  step  size  that  will  produce  accurate  results.  During  the 
initial  portion  of  the  propagation,  when  the  waveform  is  undergoing  its  most 
rapid  change  due  to  nonlinear  effects,  a  small  step  size  is  required.  After 
the  discontinuity  distance  is  reached,  however,  the  nonlinear  effects  become 
less  pronounced  and  the  step  size  can  be  increased.  In  order  to  minimize  the 
running  time  of  the  computer  program,  both  under  these  circumstances  and  on 
occasions  when  an  overly  conservative  initial  step  size  has  been  chosen, 
a  variable  step  size  feature  is  incorporated. 

This  feature  doubles  the  step  size  whenever  the  average  percentage 
change  in  the  amplitude  components  (AC^/G^I  and  | | ,  j=  1,  2,  .  .  .,  k, 
over  the  previous  step  is  below  some  arbitrary  value  e.  The  integer  k  is  the 
number  of  harmonics  printed  in  the  output  and  is  generally  less  than  the  number 
of  harmonics  retained  in  the  calculation.  There  are  two  methods  for  controlling 
the  doubling  of  the  step  size.  Either  the  internal  value  e  may  be  modified  or  the 
number  of  harmonics  printed  out  may  be  changed. 


The  uarmonics  are  printed  out  at  fixed  distance  intervals  (print  out 
distance  interval  =  specified  integer  *  initial  step  size).  When  the  step 
size  is  doubled,  it  is  unlikely  that  the  harmonics  will  be  calculated  at  positions 
coinciding  with  the  print-out  distance.  The  program  circumvents  this  problem 
by  linearly  interpolating  the  output  from  the  calculated  values. 
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D.  Renormalization 


One  of  the  standard  programing  problems  associated  with  numerical 
solutions  is  the  limited  exponent  range  of  computers.  The  Advanced  Scientific 
Computer  (ASC)  at  the  Naval  Research  Laboratory  (NRL),  for  which  this  program 
was  written,  has  an  exponent  range  of  -76  to  +76.  The  terms  most  likely  to 
exceed  this  range  are  the  exponentials  in  Eqs.  (10)  and  (11). 

To  illustrate  this  problem  the  exponentials  in  Eqs.  (10)  and  (11) 
are  examined  for  the  case  of  an  omega-squared  frequency  dependence  of  the 
absorption  coefficients  (fresh  water).  The  exponential  in  the  second  series 
on  the  right-hand  sides  of  Eqs.  (10)  and  (11)  may  then  be  written  in  the  form 

exp  [■(ak+m+am“otk)(r-r0)]  =  exp  t-2m(k+m)0t1(r-ro)] ,  (25) 

where  is  the  absorption  coefficient  of  the  fundamental.  Since  m  and  k  are 
positive  integers,  the  exponential  tends  to  zero  with  increasing  r.  This  causes 
no  problem  if  the  computer  is  told  to  set  underflow  to  zero  (the  error 
associated  with  setcing  numbers  smaller  than  10  ^  to  zero  is  negligible). 

The  exponential  in  the  first  series  on  the  right-hand  sides  of  Eqs.  (10) 
and  (11)  may  be  written  as 

exp  [-(«k+m+um-*k)(r-ro) ]  =  exp  [+2m(k-m)aj (r-rQ) ] .  (26) 

Since  m  is  always  less  than  k  in  the  first  series,  this  exponential  may  exceed 
the  upper  bound  of  the  exponent  range.  As  an  example,  when  50  harmonics  are 
retained  in  the  calculation  of  a  1 00— kHz  spherical  wave  in  fresh  water 
(ctj  =  2.38  x  10  4  N/m),  che  exponent  range  of  +76  is  exceeded  at  r  =  255  m. 

This  distance  is  totally  insufficient  to  examine  the  asymptotic  decay  of  the 
spherical  waves. 


A  simple  method  for  circumventing  this  problem  is  to  renormalize  the 
waveform  after  each  step.  In  this  procedure  the  source  position  rQ  is  changed 
after  each  step  and  set  equal  to  the  current  position  r.  This  limits  the  size 
of  the  distance  term  in  the  exponentials  to  h,  the  step  size,  and  merely  requires 
that  the  amplitudes  be  transformed  as 

\ - -Ak  [Ro/(R+h)]  exp  (-akh).  (77) 

It  has  the  additional  advantage  of  not  requiring  the  calculation  of  the  exponentials 
at  each  step  since  they  do  not  change  and  can  be  stored.  If  the  step  size  doubles, 
the  stored  values  are  simply  squared. 

IV.  DESCRIPTION  OF  COMPUTER  PROCRAM  FAW 

The  computer  program  FAW,  which  is  listed  in  Appendix  A,  is  written  in 
universal  FORTRAN  and  should  run  on  any  computer  accepting  this  language. 

However,  the  program  has  been  specifically  written  to  take  advantage  of  the 
vectorizing  capability  of  the  ASC  at  NRL. 

While  the  use  of  the  vectorizing  option  on  the  ASC  greatly  reduces  the 
running  time  of  the  program,  it  does  not  allow  the  use  of  variable  ranges  on 
nested  DO  loops.  This  constraint  results  in  FAW  running  inefficiently  on 
computers  that  do  not  have  vectorizing  capability.  Therefore,  it  may  be 
necessary  to  modify  FAW  for  use  on  a  nonvectorizing  computer. 

The  DO  loops,  which  are  modified  in  FAW  for  the  vectorizing  process, 

are  the  loops  associated  with  the  two  series  on  the  right-hand  sides  of  Eqs. 

(10)  and  (11)  (lines  243  to  256  and  320  to  333  in  FAW).  The  first  series  on 

the  right-hand  sides  of  Eqs.  (10)  and  (II)  has  an  upper  limit  of  k-1  for  the 
til 

kL  harmonic.  The  second  series  is  truncated  to  an  upper  limit  of  j-k  for  the 
kl  1  harmonic  when  j  harmonics  are  retained  in  the  calculation.  Since  both  of 
these  limits  are  a  function  of  the  harmonic  increment  being  calculated,  the 
range  of  the  inner  DO  loop  is  not  a  constant.  However,  the  ARC  runs  quicker 
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if  the  loops  are  vectorized,  and  to  do  this  the  ranges  of  the  two  inner  loops 
have  been  set  equal  to  J-l.  This  results  in  extraneous  terms  being  calculated 
and  increases  the  running  time  on  nonvectorizing  computers. 

The  extraneous  terms  do  not  contribute  to  the  calculation  and  are  all 
set  equal  to  zero  in  the  DO  loops.  This  is  accomplished  by  generating  a 
matrix  for  each  DO  loop  whose  elements  are  zero  for  the  extraneous  terms  and 
whose  nonzero  elements  are  the  exponential  terms  in  Eqs.  (10)  and  (11).  As 
explained  in  the  section  on  renormalization,  this  matrix  need  only  be  calculated 
once  and  is  then  used  at  each  step. 

The  remainder  of  this  section  describes  FAW.  Included  are  a  listing 
of  the  significant  FORTRAN  variable  names  and  descriptions  of  the  major 
computation  blocks  in  FAW  followed  by  a  discussion  of  the  parameter  input, 
printed  output,  number  of  harmonics  to  retain  in  calculation,  and  initial 
step  size. 

A.  Significant  FORTRAN  Variable  Names 

The  significant  FORTRAN  variable  names  in  FAW  are  as  follows: 

A  Geometrical  spreading  factor 

=  0  (plane  waves), 

=  1/2  (cylindrical  waves) , 

=  1  (spherical  waves). 

tVi 

ALPHA  Vector  whose  1  element  ALPHA  (I)  is  the  absorption 

coefficient  of  the  1  harmonic. 

B  The  constant  b  =  (5/c  ^  in  Eqs.  (10)  and  (11). 


BETA 


Coefficient  ^f  nonlinearity 


Small  signal  sound  speed  cq. 


Vector  containing  the  differential  change  of  the  amplitude 
coefficients  G.  The  element  DG(I)  contains  the  differential 
change  hG^(N)  at  the  N  step  as  given  hy  Eq.  (21). 


Vector  whose  IL  element  DG1(I)  contains  the  differential 
*  t*  h 

change  hCj(f’+l)  at  the  N+1  step  as  given  by  Fq.  (??)• 


Vector  containing  the  differential  change  of  the  amplitude 
coefficients  H.  The  element  Dll ( I )  contains  the  differential 
change  hll  (N)  at  the  step  as  given  by  Eq.  (22). 


Vector  whose  I  element  DU 1(1)  contains  the  differential 
change  hli'(N+l)  at  the  K+l  step  as  given  by  Eq.  (24). 

Current  step  size. 

Initial  step  size. 

Normalization  constant.  The  output  is  normalized  to  the 
constant  E, which  is  an  input  parameter  and  generally  set 
equal  to  the  initial  amplitude  of  the  fundamental. 


Factor  from  Eqr.(10)  and  (11)  equal  to  (rQ/r)  . 


Frequency  of  fundamental. 


Vector  whose  Ic  element  contains  the  amplitude  C(I)  of 

t  h 

the  sine  component  of  the  I  harmonic. 


Vector  whose  I  element  G1I(J)  =  exp  [-ALPHA  (I)  *  DX]. 
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Array  whose  elements  are  all  zero  that  is  used  as  a  buffer 
for  the  working  array  02.  In  the  DO  loops  which  calculate 
DC,  DG1,  DH,  and  DH1  the  ranges  of  the  loops  have  been 
written  to  take  advantage  of  vectorization.  This  results 
in  negative  indices  for  some  of  the  02  elements  in  the 
calculation.  The  GX  array  is  placed  before  the  G2  array  to 
prevent  incorrect  results. 

til 

Array  whose  I  element  first  contains  G(I)  and  in  later 
calculations  contains  G(l)  +  DG(I)  as  required  by  the 
Runge-Kutta  method. 

Same  as  G  for  the  cosine  elements. 

Same  as  GX  but  placed  in  front  of  112. 

Same  as  G2  for  the  cosine  elements. 

Input  parameter  equal  to  either  zero,  if  only  the  absorption 

coefficient  for  the  fundamental  is  input  and  an  omega- 

squared  dependence  for  the  harmonic  absorption  coefficients 
is  used,  or  one  if  the  absorption  coefficient  is  input  for 
each  harmonic  retained  in  the  calculations. 

Print-out  interval  =  IP  *  DX. 

Number  of  harmonics  retained  in  the  calculations. 

Number  of  input  G  coefficients. 

Number  of  input  H  coefficients. 


1 A 
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Number  of  externally  supplied  absorption  coefficients 
used  to  modify  the  omega-squared  dependent  coefficients. 

rrly  if  IAF  =  0  and  the  omega-squared  dependence  is  to 
be  modified. 

KI  Number  of  harmonics  printed  to  output. 

R  Spatial  coordinate. 

RMAX  Maximum  distance  that  propagation  is  to  be  calculated. 

RN  Source  position.  Modified  after  each  step.  Set  equal  to 

the  present  position  R  as  described  in  section  on  renormal- 
izat ion. 

RO  Initial  source  position. 

XI  Doubly  dimensioned  array  whose  elements  XI  (M,  N)  are  unity 

if  XI 1  (M,  K)  is  nonzero  and  zero  otherwise. 

XI I  Array  whose  elements  XI 1  (M,  N)  =  exp  ((ALPHA(M)  - 
ALPHA(N)  -  ALPHA  (M-N))  *  DX)  for  M  =  ?  to  J  and  N  =  1 
to  M.  All  other  elements  are  zero. 

XZ  Doubly  dimensioned  array  whose  elements  X2  (M,  N)  are  unity 

if  X22  (M,  K)  is  nonzero  and  zero  otherwise. 

X22  Array  whose  elements  X22  (M,  N)  =  exp  ((  ALPHA(M)  - 

ALPHA(N)  -  ALPHA  (M+N))  *  DX)  for  M  =  1  to  J-l  and 
K  =  1  to  J-M.  All  other  elements  are  zero. 
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B.  Major  Computation  Blocks  in  FAW 


Descriptions  of  the  major  computation  blocks  in  FAW  are  as  follows: 


Computation  Block 


Initialize  all  arrays  to  zero. 


Line  Number 
From  To 

29  44 


Read  first  fourteen  data  cards  and 
print  data  to  output. 


51  138 


Obtain  absorption  coefficients. 


Read  input  waveform  and  print  waveform 
to  output. 


142  179 


183  199 


Calculate  matrix  elements  of  Xl, 
XI 1,  X2,  and  X22  and  elements  of 
vector  GH. 


204  221 


Calculate  DG  and  DH,  first  derivative 
of  Runge-Kutta  method. 


243  256 


Calculate  new  amplitudes  for  second 
derivative. 


260  265 


Find  last  five  non-zero  harmonic 
amplitudes  and  modify  them,  if 
necessary,  to  insure  that  they  form 
a  non-increasing  sequence. 

Calculate  DG1  and  DH1,  second 
derivative  of  Runge-Kutta  method. 


289  308 


320  333 


IH*-  '^sIB^Siaisawww 
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Computation  Block 

Line 

Number 

From 

To 

Calculate  new  amplitudes  for  current 

position. 

337 

340 

Find  last  five  non-zero  harmonic 

amplitudes  and  modify  them,  if  necessary, 

to  insure  that  they  form  a  non-increasing 

sequence. 

344 

382 

Check  position  for  output. 

386 

386 

Interpolate  output  if  output  position 

does  not  coincide  with  position  of 

calculated  amplitudes. 

390 

404 

Standard  output  (no  interpolation 

necessary). 

405 

408 

Print  output. 

409 

418 

Check  step  size  (if  incremental 

change  is  small,  then  double  step 

size). 

422 

431 

Modify  waveform  for  renormalization. 

435 

438 

Double  step  size. 

443 

454 

Next  step. 

456 
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C.  Parameter  Input 


The  input  consists  of  a  series  of  data  cards: 

Data  Card  I:  Format  D17.10  -  FREQ,  the  frequency  of  the  fundamental. 

Data  Card  2:  Format  D17.10  -  C,  the  small  signal  sound  speed  in 
meters  per  second. 

Data  Card  3:  Format  D17.10  -  BETA,  the  coefficient  of  nonlinearity. 

Data  Card  4:  Format  D17.10  -  E,  the  normalization  constant. 

Data  Card  5:  Format  D17.10  -  DX1,  the  initial  step  size  in  meters. 

Data  Card  6:  Format  D17. 10  -  RMAX,  the  maximum  propagation  distance 
in  meters. 

Data  Card  7:  Format  D17.10  -  RO,  the  source  size,  in  meters,  for 

cylindrical  and  spherical  waves.  If  a  plane  wave  is  being 
calculated,  RO  is  not  used  but  a  value  must  be  entered. 

Data  Card  8:  Format  F5. 2  -  A,  the  geometrical  spreading  factor. 

A  =  0  (plane  waves) 

=  1/2  (cylindrical  waves) 

=  1  (spherical  waves). 

Data  Card  9:  Format  14  -  IAF,  the  absorption  flag.  IAF  =  0  if  the 

omega-squared  dependence  for  the  absorption  coefficients  is 
being  used.  IAF  =  1  if  all  absorption  coefficients  are 
entered  on  data  cards. 

Data  Card  10:  Format  14  -  .J,  the  number  of  harmonics  retained  in 
the  calculation. 


Data  Card  11:  Format  14  -  K,  the  number  of  initial  G  coefficients  entered. 

Data  Card  12:  Format  14  -  KI,  the  number  of  initial  H  coefficients  entered. 

Data  Card  13:  Format  14  -  NI,  the  number  of  harmonics  printed  to  output. 

Data  Card  14:  Format  18  -  IP,  the  integer  multiplicative  factor  of 

DX1  which  gives  the  print-out  interval. 

If  IAF  =  0 

Data  Card  15:  Format  D17.10  -  ALPHA  (1),  the  absorption  coefficient 
of  the  fundamental  in  nepers  per  meter. 

Data  Card  16:  Format  14  -  L,  the  number  of  harmonics  being  modified 
from  the  omega— squared  dependence. 

Next  L  Cards:  Format  14,  D17. 10  -  These  cards  contain  an  integer,  right 

justified  in  the  first  four  spaces  on  the  card,  specifying  the 
number  of  the  harmonic  followed  by  the  absorption  coefficient 
for  that  harmonic. 

If  IAF  =  1 

Data  Card  15  to  14+J:  Format  D17. 10  -  These  cards  contain  the  J  absorption 
coefficients  In  order  in  nepers  per  meter. 

Next  K  Cards:  Format  14,  D17. 10  -  These  cards  contain  a  right  justified 

integer,  in  the  first  four  spaces,  specifying  the  number  of 
the  harmonic  followed  by  the  G  amplitude  coefficient  of  that 
harmonic  in  meters  per  second. 

Next  ”T  Carls:  Format  14,  D17. 10  -  These  cards  contain  a  right  justified 

integer,  in  the  first  four  spaces,  specifying  the  number  of 
the  harmonic  followed  by  the  H  amplitude  coefficient  of 
that  harmonic  in  meters  per  second. 
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D.  Printed  Output 


A  sample  output  from  FAW  is  shown  in  Appendix  R.  The  first  portion  of 
the  output  contains  a  listing  of  all  input  parameters.  This  allows  the  input 
parameters  to  be  checked  for  errors  and  is  useful  for  future  reference.  The 
remaining  portion  contains  the  calculated  output. 

The  distance  in  meters  is  printed  at  each  output  interval  followed  by 
two  columns  of  numbers.  The  columns  are  labeled  and  contain  the  amplitudes 
of  the  sine  and  cosine  components  divided  by  the  normalization  constant  E. 

The  number  of  terms  output  at  each  interval  is  NI,  which  is  a  user-specified 
parameter.  In  addition,  the  step  size  AX  in  meters,  is  printed  each  time  the 
step  size  is  doubled.  This  gives  the  user  useful  information  on  the  effects 
of  varying  the  step  size  doubling  parameter  r  and  the  number  of  harmonics 
printed  out  NI. 

E.  Number  of  Harmonics 

A  sufficient  number  of  harmonics  must  be  retained  in  the  calculation  to 
insure  a  negligibly  small  error  in  the  highest  harmonic  of  interest.  In  order 
to  obtain  some  measure  of  the  required  number  of  harmonics,  the  algorithm 
was  used  to  calculate  the  harmonic  content  of  an  initially  pui e  sinusoidal 
plane  wave  with  a  variety  of  harmonics  being  retained  in  the  calculation.  The 
frequency  of  the  fundamental  was  2.5  MHz,  and  the  initial  pressure  amplitude 
was  3  atmospheres,  which  gives  a  discontinuity  distance  of  21  cm.  Table  I 
lists  the  percentage  deviations  of  the  resulting  amplitudes  of  the  first  five 
harmonics,  at  the  discontinuity  distance,  from  the  values  obtained  when  40 
harmonics  were  retained.  As  is  obvious  from  the  table,  the  required  number 
of  harmonics  to  retain  depends  on  the  harmonic  of  interest  and  the  allowable 
error.  If  the  fundamental  is  the  only  harmonic  of  interest,  one  need  retain 
no  more  than  seven  harmonics  in  the  calculation.  On  the  other  hand,  an 
accurate  value  for  the  fifth  harmonic  may  require  twenty  or  more  harmonics 
to  be  retained  in  the  calculation. 


Table 

1  -  !  c roinaye 
for  various 

leviation  in 
ntimLurs  of 

the  first  five  hartt’ii 
retained  harmonics  M 

ics 

II/.KIOHIC^w 

7 

8 

9 

10 

15 

20 

25 

30 

1st 

b.03 

0.<»2 

0.02 

0.02 

0.004 

0.0003 

0.00002 

0.00000 

2nd 

0.20 

0.20 

0. 10 

0.10 

0.020 

0.0020 

0.00030 

0.00003 

3rd 

0.80 

0.60 

0.50 

0.50 

0.080 

0.0090 

0.00100 

0.00010 

Hill 

2.40 

1.80 

1.60 

1.50 

0.200 

0.0200 

0.00100 

0.00040 

Uh 

4.40 

3.80 

3.50 

3.40 

0.500 

0.0500 

0.00700 

0.00080 

F.  Step  Size 

With  the  ability  of  the  program  to  double  the  step  size  after  each  step, 
it  is  best  to  choose  the  initial  step  size  conservatively  and  let  the  program 
find  the  best  value.  In  order  to  determine  a  conservative  initial  valtie  for 
the  step  size,  the  choice  of  step  size  was  investigated  for  the  same  2. S-MHz 
plane-wave  case  used  in  the  previous  section.  The  step  size  was  not  allowed 
to  double  and  various  step  sizes  from  1/200  to  1/10  of  the  discontinuity 
distance  were  used.  Table  II  lists  the  pe  centage  deviations  of  the  first 
five  harmonic  amplitudes  from  the  values  obtained  when  the  step  size  was  1/200 
of  the  discontinuity  distance.  The  amplitudes  were  those  at  the  discontinuity 
distance,  and  forty  harmonics  were  retained  in  the  calculations.  The  table 
indicates  that  a  step  size  of  1/10  of  the  discontinuity  distance  will  yield 
results  that  are  accurate  to  within  the  normal  experimental  error. 


Table  II  -  Percentage  deviation  in  the  first  five  harmonics 
for  various  step  sizes 

(o  =  discontinuity  distance) 


-^STEP  SIZE 

IARMONIcT' 

1/100  a 

1/5C  o 

1/20  o 

1/10  o 

1st 

0.0005 

0.0030 

0.0200 

0.0600 

2nd 

0.0001 

0.0050 

0.0300 

0.1000 

3rd 

0.0020 

0.0100 

0.0900 

0.4000 

4  th 

0.0070 

0.0400 

0.2500 

1.1000 

5th 

0.0100 

0.0700 

0.5000 

2.0000 

V.  COMPARISON  WITH  PHENOMENOLOGICAL  MODEL 


As  a  test  of  their  validity,  Eqs.  (10)  and  (11)  were  used  to  compute  the 
harmonic  content  of  an  initially  pure  sinusoidal  300-Hz  plane  wave.  An  omega- 
squared  frequency  dependence  of  the  linear  absorption  terns  was  assumed  with 
the  absorption  of  the  fundamental  being  1.12  x  10~fiNp/m.  The  initial  pressure 
was  1  atmosphere.  This  problem  was  also  solved  using  a  computer  algorithm 
[6]  based  on  the  phenomenological  model  of  Fox  and  Wallace  [7].  Figure  (1) 
illustrates  the  agreement  between  the  results  obtained  using  the  algorithm 
presented  in  this  paper,  shown  as  solid  curves,  and  the  results  obtained  using 
the  phenomenological  model,  shown  as  dots.  Although  only  the  first  four 
harmonics  are  illustrated  in  this  figure,  the  agreement  was  equally  as  good 
for  higher  harmonics. 
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3RD  HARMONIC 


DISTANCE  (km) 


Calculated  harmonic  behavior  of  an  initially  pure^. 
sinusoidal  plane  wave  (v  =  300  Hz,  *  =  1.12  x  10  N/m. 

and  P  =  1  atm.)  by  -  frequency  domain  algorithm, 

•  phenomenological  model. 
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COMPUTER  PROGRAM  FAW 
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THIS  PROGRAM  STEPWISE  CALCULATES  THE  HARMONIC  CONTENT 
OF  A  FINITE  AMPLITUDE  WAVE  AS  A  FUNCTION  OF  POSITION. 
THE  INITIAL  WAVEFORM  MAY  HAVE  ANY  HARMONIC  CONTENT 
AND  ARBITRARY  PHASE.  THE  ABSORPTION  COEFFICIENTS  MAY 
BE  EXTERNALLY  SUPPLIED  OR  AN  INTERNAL  OMEGA  SQUARED 
ALGORITHM  IS  SUPPLIED.  THE  INTERNAL  ALGORITHM  ALLOWS 
MODIFICATION  OF  ANY  OF  THE  OMEGA  SQUARED  COEFFICIENTS. 
THE  NUMERICAL  INTEGRATION  USES  THE  RUNGE-KUTTA  METHOD 
AND  THE  STEP  SIZE  WILL  AUTOMATICALLY  INCREASE  WHEN  THE 
CHANGE  OVER  AN  INTERVAL  IS  SMALL. 


DIMENSION  G  ( 50 ) »  GX ( 50 ) >62(100) »H(50) » HX ( 50 ) » H2  ( 100) 
DIMENSION  XI  (50 >50) > XI 1(50 >50) >X2(50>50) 

DIMENSION  X22(50»50) >  XL ( 50 ) >  XL1 (50) >  GH ( 50 ) 

DIMENSION  ALPHA (50) > DG(50> > DG 1 ( 50 ) > DH ( 50 ) >DH1(50) 
DIMENSION  RG  ( 5 ) >  KH ( 5 ) 

DOUBLE  PRECISION  G>GX?G2>H>HX>H2>X1>X11>X2>X22>XL>XL1 
DOUBLE  PRECISION  GH > ALPHA > DG > DG 1 > DH f DH1 > R > C > FREQ > BETA 
DOUBLE  PRECISION  E  >  B  >  X>  FI >  DX  >  RMAX  >  RN  >  Z  >  Z1 >  DXI >PD 
DOUBLE  PRECISION  YrRO 

INITIALIZE  REGISTERS  TO  ZERO 

DO  110  IB=1 > 100 
G ( I B ) =0 . DO 
GX ( I B ) =0 . DO 
G2 ( I B ) =0 . DO 
H  ( I B )  =  0  .  D  0 
HX ( I B ) =0 . DO 
H2(  I B ) - 0 , D 0 
CONTINUE 
DO  ISO  IC  =  1 >  50 
DO  120  J  Ii  - 1  r  50 

XI  (IC»  110=0.  DO 

XI I  ( IC»  ID>=0.)i0 
X2 ( IC » I D ) =0 . DO 
X22 ( IC » f  D ) =0 , DO 
CONTINUE 

continue: 

INPUT  DATA 


FREQ-FUNDAMENTAL  FREQUENCY 


READ ( 5  >  1 40 ) FREQ 


25 


52 

53 

140 

54 

150 

55 

C 

56 

C 

57 

58 

59 

C 

60 

160 

61 

C 

62 

C 

63 

64 

65 

C 

66 

170 

67 

C 

68 

C 

69 

70 

71 

C 

72 

180 

73 

C 

74 

C 

75 

76 

77 

C 

78 

190 

79 

C 

80 

C 

81 

82 

83 

C 

84 

200 

85 

c 

86 

c 

87 

88 
89 

C 

90 

210 

91 

c 

92 

c 

93 

c 

94 

c 

95 

i; 

96 

97 

c 

98 

99 

220 

100 

230 

101 

c 

102 

c 

FORMAT ( D17 . 10  ) 

PRINT  150  , FREQ 

FORMAT ( 5X, 10HFREQUENCY® , D17 . 10) 

C=SMALL  SIGNAL  SOUND  SPEED 

READ<5> 140)C 
PRINT  160 ,  C 

FORMAT ( 5X , 12HS0UND  SPEED® »D 17 .10) 

BETA=COEFFICIENT  OF  NONLINEARITY  ( 1  +  B/2A  ) 

READ(5» 140)BETA 

PRINT  170, BETA 

FORMAT ( 5X , 5HBET  A=»D17.10) 

E=NORMALIZATION  CONSTANT  (OUTPUT  IN  DB  RE ( E )  ) 

READ ( 5 , 1 40 ) E 
PRINT  180, E 

FORMAT (5X,23HN0RMALIZATI ON  CONSTANT®, D 17. 10) 

DXI=INITIAL  STEP  SIZE 

READ  <  5 , 1 40 ) DXI 
PRINT  190, DXI 

FORMAT <5X»18H INITIAL  STEP  SIZE® , D17 . 10 ) 

RMAX=MAXIMUM  PROPAGATION  DISTANCE 

READ (5,140) RMAX 
PRINT  200, RMAX 

FORMAT (5X, 17HMAXIMUM  DISTANCE® , D1 7 . 10 ) 

RO=SOURCE  SIZE( INITIAL  POSITION  FOR  GEOMETRICAL  SPREADING) 

READ ( 5 , 1 40 ) RO 
PRINT  210, RO 

FORMAT ( 5X, 12HS0URCE  SIZE® , D17 . 10 > 

A=SPREADING  FACTOR 
=0  (PLANE  WAVES) 

=1/2  (CYLINDRICAL  WAVES ) 

=1  (SPHERICAL  WAVES) 


READ ( 5 , 220 ) A 
FORM AT (F5. 2) 

PRINT  230,  A 

FORMAT ( 5X, 17HSPREADING  F ACTOR® , F5 . 2 ) 
IAF-  FLAG (ABSORPTION  COEFFICIENTS  ) 


103 

104 

105 

106 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 


=0  (OMEGA  SQUARED  DEPENDENCE) 

=1  (EXTERNALLY  SUPPLIED) 

READ ( 5  , 240 ) I AF 
FORMAT ( 14) 

PRINT  250 ,  I AF 

FORMAT (5X,16HABS0RPTI0N  FLAG*, 14) 

J=NUMBER  OF  HARMONICS  RETAINED  IN  CALCULATION 

READ ( 5 »  240 ) J 
PRINT  260,  J 

FORMAT ( 5X  »  20HNUMBER  OF  HARMONICS* » 14 ) 

K=NUMBER  OF  INITIAL  G  COEFFICIENTS 
K I “NUMBER  OF  INITIAL  H  COEFFICIENTS 

READ ( 5 » 240 ) K 
READ ( 5  ,  240 ) KI 
PRINT  270, K 

FORMAT (5X»23HINITIAL  G  COEFFICIENTS*, 14 ) 

PRINT  280, K I 

F0RMAT(5X»23H INITIAL  H  COEFFICIENTS*, 14 ) 

NI=NUMBER  OF  HARMONICS  PRINTED  OUT 

READ ( 5 , 240 ) NI 
PRINT  290, NI 

FORMAT ( 5X , 27HNUMBER  OF  HARMONICS  OUTPUT*, 14) 

IP-PRINT  OUT  INTERVAL* I P*DX 

READ ( 5 , 300 ) IP 
FORMAT( 18) 

PRINT  310, IP 

FORMAT ( 5X , 19HPRINT  OUT  INTERVAL* , IS ) 

OBTAIN  ABSORPTION  COEFFICIENTS 

IF(IAF)1110,320,400 

OMEGA  SQUARED  DEPENDENCE 

READ ( 5 , 140 ) ALPHA ( 1 ) 

PRINT  330 , ALPHA (  1 ) 

F0RMAT(5X,9HALPHA( 1 )  =  *  D17 . 1 0 ) 

DO  340  I E* 2  ,  .1 

ALPHA ( IE )=IE*IE*ALPHA(1 ) 

CONTINUE 

L*NUMBER  OF  ABSORPTION  COEFFICIENTS  BEING  MODIFIED 


FROM  OMEGA  SQUARED  DEPENDENCE 


154  C 

155  C 

156  READ  <  5 »  240  )  L 

157  PRINT  350, L 

158  350  FORMA  r ( 5X , 43HNUMDER  OF  MODIFIED  ABSORPTION  COEFFICIENTS* » 14 ) 

159  IF(L) 1110,440,360 

160  360  PRINT  370 

161  370  FORMAT  <  5X , 33HM0DIFIED  ABSORPTION  COEFFICIENTS.) 

162  DO  390  IG=1 ,L 

163  READ ( 5 , 380 ) N , X 

164  380  FORMAT! 14, D17. 10) 

165  PRINT  380, N,X 

166  ALPHA(N)=X 

167  390  CONTINUE 

168  GO  TO  440 

169  C 

170  C  INPUT  ABSORPTION  COEFFICIENTS 

171  C 

172  400  PRINT  410 

173  410  FORMAT  <5X»23HABS0RPTIQN  COEFFICIENTS) 

174  DO  430  IH= 1 » J 

175  RE AD (5, 140)X 

176  ALPHA ( IH ) =X 

177  PRINT  420»IH»ALPHA<IH) 

178  420  FORMAT ( 1X»6HALPHA! , I3,2H)=»D17.10) 

179  430  CONTINUE 

180  C 

181  C  READ  INPUT  WAVEFORM 

182  C 

183  440  CONTINUE 

184  IF!I\)  1110,490,450 

185  450  PRINT  460 

186  460  FORMAT (5X,14HINFUT  WAVEFORM) 

187  DO  480  I J=1 ,  K 

188  READ ( 5 , 380 ) N  ,  X 

189  G ( N ) =  X 

190  PRINT  470 , N , G ( N ) 

191  470  FORMAT! 1X»2HG! , 13, 2H)=,D17. 10) 

192  480  CONTINUE 

193  490  IF(KI) 1110,530,500 

l<?4  500  DO  520  IK=1,K1 

195  READ ( 5 , 380 ) N , X 

196  H(N)=X 

197  PRINT  510»N,HlN> 

198  510  FORMAT! 1X,2HH( , I3,2H)=,D17. 10) 

199  520  CONTINUE 

200  530  B=3.1415926536*BETA*FREQ/!C*C) 

201  C 

202  C  CALCULATE  MATRIX  ELEMENTS 

203  C 

204  DO  550  IL  =  2 » J 


28 


1 

205 

IL1=IL-1 

1 

206 

DO  540  IM=1 » IL1 

1 

207 

Xl.t  (  I.L>  IM)  =DEXP <  ( ALPHA <  IL) -  ALPHA <  IM > -  ALPHA ( IL-IM )  )  *DXI ) 

ii 

i 

208 

XI ( IL » IH )  =  1 . DO 

i 

H 

209  540 

CONTINUE 

R 

1 

210  550 

CONTINUE 

se 

i 

211 

IFX= J-l 

| 

212 

DO  570  IN=1 f IPX 

£ 

213 

IN1=J-IN 

i :  * 

214 

DO  560  IP=1 » INI 

215 

X22 ( IN» IP > =DEXP < (ALPHA (IN) -ALPHA (IP) -ALPHA (IN+IP) )*DXI ) 

216 

X2 ( IN  f IP)  =  1 . DO 

s 

217  560 

CONTINUE 

; 

218  570 

CONTINUE 

i 

219 

DO  580  IQ=lfJ 

220 

GH ( IQ ) =DEXP ( -ALPHA  < IQ ) #DXI ) 

221  580 

222  C 

CONTINUE 

$ 

223  C 

224  C 

SET  COUNTER 

1 

225 

RN=RO 

» 

f 

226 

R=RN 

| 

227 

DX=DXI 

228 

PD=IP*DX 

i 

229 

Fl=l .DO 

1 

230  590 

CB=DX*B 

i 

231  C 

i 

232  C 

233  C 

FILL  ARRAYS 

i 

234 

DO  600  I R= 1 »  J 

I 

235 

G2(IR)=G(IR) 

236 

H2 ( IR ) =H ( IR ) 

c 

1 

237 

DG ( IR ) =0 . DO  ' 

i 

238 

DH ( IR ) =0 . DO 

239  600 

240  C 

CONTINUE 

! 

241  C 

242  C 

ENTER  LOOP  FOR  CALCULATING  FIRST  DERIVATIVE 

243 

DO  620  IS=2»J 

j 

244 

DO  610  IT=1 » IFX 

245 

CX=CB*IT*X1(IS»IT) 

j 

246 

DG(IS)=CX*(G2<IS-IT)*G2<IT)-H2(IS-IT)*H2(IT))+DG<IS) 

I 

247 

DH( IS)=CX*(H2( IS-IT)*G2( IT)+G2( IS-IT)*H2( IT) )+DH( IS) 

248  610 

CONTINUE 

f 

249  620 

CONTINUE 

1  * 

250 

DO  640  IU=1»J 

1 

1 

251 

DO  630  IV=1,IFX 

i 

252 

CX  =  CBKIll*X2(IU»IV)  , 

|  * 

253 

DG<  IU)=DG(  IU)-CX»(G2aU+IV)*G2(IV)+H2(IU+IV)*H2<  IV) ) 

i 

254 

DH< IU)=DH( IU)+CX*(G2< IU+I V ) »H2< I V ) -H2( IU+I V ) *G2 ( I V ) ) 

I 

1 

255  630 

CONTINUE 

i 


i 

i 

i 

i 

K 

I 

s 


29 


256 

640 

257 

C 

258 

C 

259 

C 

260 

261 

262 

263 

264 

265 

650 

266 

C 

267 

C 

268 

C 

269 

270 

271 

272 

655 

273 

274 

275 

660 

276 

277 

670 

278 

279 

680 

280 

690 

281 

282 

700 

283 

284 

710 

285 

C 

286 

C 

287 

C 

288 

C 

289 

290 

291 

29° 

712 

293 

294 

295 

714 

296 

297 

298 

720 

299 

300 

301 

302 

790 

303 

304 

305 

724 

306 

CONTINUE 

CALCULATE  NEW  AMPLITUDES  FOR  SECOND  DERIVATIVE 

DO  650  IW=1»J 
G2(IW)=G(IW)+DG<IW> 

H2(IW)=H(IW>+DH<IW> 

DG 1 ( I W ) =0 . DO 
DH1 ( I W ) =0 . DO 
CONTINUE 

FIND  LAST  FIVE  NON-ZERO  HARMONICS 

DO  655  IW1  =  1 >  5 
KG ( IU1 ) =0 
KH ( IU1 ) =0 
CONTINUE 
DO  710  IX=1,J 
IF(G2(IX>  >660.680.660 
DO  670  I Y=1 *  4 
KG ( I Y ) =KG ( I Y  + 1 ) 

CONTINUE 
KG ( 5 )  =  I X 

IF(H2( IX) ) 690 .710.690 
DO  700  IZ=1 » 4 
KH(IZ)=KH(IZ+1) 

CONTINUE 
KH ( 5 ) = I X 
CONTINUE 

INSURE  THAT  THE  LAST  FIVE  HARMONICS  ARE  NOT 
PROGRESSIVELY  LARGER 

DO  720  JA= l ?  4 
JB -KG ( JA ) 

IF (JB)712»720»712 
Z=DABS(G2(JB> ) 

JC  =  KG ( JA+1 ) 

IF(JC)714.720.714 
Z1  =  DABS ( G2( JC  )  ) 

IF ( Z1 . LT . Z ) GO  TO  720 

G2(  JC) =0.95*02 <JC)*Z/Z1 

CONTINUE 

DO  730  JD=1 »  4 

JE=KH ( JD ) 

IF< JE>722,730.722 
Z=DABS(H2< JE) ) 

JF =KH ( JD+1 ) 

IF(JF >724.730.724 
Z1  =  DABS(H2( JF) ) 

IF ( Z1 «  LT . Z ) GO  TO  730 


3C 


307 

308  730 

309 

310 

311  740 

312 

313 

314 

315  750 

316 

317  C 

318  C 

319  C 

320  760 

321 

322 

323 

324 

325  770 

326  780 

327 

328 

329 

330 

331 

332  790 

333  800 

334  C 

335  C 

336  C 

337 

338 

339 

340  810 

341  C 

342  C 

343  C 

344 

345 


H2< JF)=0.95*H2< JF)*Z/Z1 

CONTINUE 

R=R+DX 

<IF  <  A  )  1 1 10  *  760  f  740 
IF ( A . LT .  1 . DO ) GO  TO  750 
F1=RN/R 
CB=DX*B*F1 
GO  TO  760 
F1=<RN/R)¥*0.S 
CB=DX*B*F1 

ENTER  LOOP  FOR  CALCULATING  SECOND  DERIVATIVE 

DO  780  JG=2 »  J 

DO  770  JH=1 » IFX 
CX=CB#JH#X11<JG»JH) 

DG1( JG)=CX*<G2< JG-JH)*G2< JH)-H2( JG-JH)*H2< JH) >+DGl< JG> 
DH1 ( JG ) =CX# (H2( JG- JH ) #G2( JH>  +G2<  JG-JH ) *H2 ( JH ) ) +DH1 ( JG ) 
CONTINUE 
CONTINUE 
DO  800  JI=lfJ 

DO  790  JJ=1 » IFX 
CX=CB*JI*X22< JI, JJ> 

DG1( JI )=DG1( JI)-CX*(G2< JI+JJ>*G2( JJ)+H2C JI+JJ)*H2< JJ) ) 
DH1 <  JI ) =DH1 ( JI )+CX*(G2<  JI  + JJ)*H2( JJ) -H2 ( JI+JJ )#G2 ( JJ ) > 
CONTINUE 
CONTINUE 

CALCULATE  NEW  AMPLITUDES 
DO  810  JK=1 / J 

G ( JK  >  =G  <  JK ) +0 . 5* ( DG ( JK ) +DG 1 ( JK  > ) 

H  (  JK )  =H  ( JK  )  +0 . 5*  ( DH  ( JR )  +DH1  (  JK  )  ) 

CONTINUE 

FIND  LAST  FIVE  NON-ZERO  HARMONICS 

DO  320  JL= 1 >  5 
KG( JL)=0 


346  KH ( JL ) =0 

347  320  CONTINUE 

348  DO  8S0  JM=1 »  J 


349  IF(G(JM) )830»850f330 

350  830  DO  B40  JN=1»4 

351  KG(JN)=KG( JM+1) 

352  840  CONTINUE 


353 

KG ( 5 ) = JM 

354 

850 

IF ( H ( JM ) ) 860 1 880  *  860 

355 

360 

DO  870  JG=1 f  4 

356 

KH ( JO ) =KH ( J0+ l ) 

357 

870 

CONTINUE 

?! 


358  KH ( 5 ) = JM 

359  830  CONTINUE 

360  C 

361  C  INSURE  THAT  THE  LAST  FIVE  HARMONICS  ARE  NOT  PROGRESSIVELY  LARGER 

362  C 

363  DO  890  JP=1»4 

364  JQ=KG ( JP) 

365  IF< JQ)882.890.882 

366  382  Z=DABS(G<  JO)  ) 

367  JR=KG(  JP  +  1 ) 

368  IF (JR >884 .390 .884 

369  884  Z1=DABS<G< JR) ) 

370  IF ( Z 1 . LT . Z) 60  TO  390 

371  G<  JR)=0.95*G<  JR)*Z/Z1 

372  890  CONTINUE 

373  DO  900  JS=  1  *  4 

374  JT=KH ( JS ) 

375  IF(JT)892»900»892 

376  892  Z=DABS(H< JT> > 

377  JU=KH( JS  +  1 ) 

378  IF< JUJ894.900. 394 

379  894  Z 1  =  DABS  <  H ( J U ) ) 

380  IFIZI.LT.ZIGO  TO  900 

381  H<  JIJ)=0.95*H<  JU)*Z/Z1 

382  900  CONTINUE 

383  C 

384  C  CHECK  FOR  OUTPUT  DISTANCE 

385  C 

386  905  IF( (R-RO)-PD) 1010.950.910 

387  C 

388  C  INTERPOLATE  OUTPUT 

389  C 

390  910  Y=DX-< (R-RO)-PD) 

391  DO  «20  JV=1»NI 

392  G2< JV)=G2< JV)-DG< JV)+0.5*(DG( JV1+DG1 < JV>>*Y/DX 

392  H2( JV)=H2( JV)-DH( JV)+0.5*(DH( JVI+DHK JV) )*Y/DX 

394  920  CONTINUE 

395  X=F 1 

396  IF(A-0. 5)940. 925. 930 

397  925  X=(RN/(RN+Y) )**0.5 

398  GO  10  940 

399  930  X=RN/ ( RN+ T ) 

400  940  DO  945  JV1-1.NI 

401  XL  (  JV1 )  -X*G2<  JV1  >  *DEXP< -ALF'HAI  JV1  )*Y)/E 

402  XL  1 ( JV 1 ) =X*H2 ( JV 1 ) *DEXF ( -ALPHA ( JV1 >*Y )/E 

403  945  CONTINUE 

404  GO  TO  965 

405  950  DO  960  JX-l.Nt 

406  XL ( JX ) - F 1 *G <  JX ) *0H ( JX ) /E 

407  XL1 < JX)=Fl*H< JX)tGH( JX)/E 

408  960  CONTINUE 

\ 


409 

965 

PRINT  970 » PD 

410 

970 

FORMAT <5X,9HDISTANCE=,D17. 10) 

411 

PRINT  980 

412 

980 

FORMAT (5X,3HSIN»20X,3HC0S) 

413 

DO  1000  JY=1,NI 

414 

PRINT  990 , JY , XL ( JY  ) , XL1 ( JY ) 

415 

990 

F0RMAT<1X»I4»4X*D17.10»5X,D17.10) 

416 

1000 

CONTINUE 

417 

PD=PD+IP*DXI 

413 

GO  TO  905 

419 

C 

420 

c 

CHECK  STEP  SIZE 

421 

c 

422 

1010 

N=0 

423 

Z  =  0 .  DO 

424 

DO  1050  J 1  =  1 »  N I 

425 

IF<G<J1) >1020,1030,1020 

426 

1020 

Z=Z+DABS<0.5*< (D6< J1 ) +DG1 (J1))/G(J1>) ) 

427 

N=N+1 

428 

1030 

IF(H(J1) ) 1040, 1050, 1040 

429 

1040 

Z=Z+DABS (0*5#(DH(J1) +DH1 (J1))/H(J1)) 

430 

N=N  +  1 

431 

1050 

CONTINUE 

432 

C 

433 

C 

MODIFY  WAVEFORM 

434 

C 

435 

DO  1060  J2=1.J 

436 

G( J2)=G< J2)*F1*GH< J2) 

437 

H( J2)=H( J2)*F1*GH< J2) 

438 

1060 

CONTINUE 

439 

IF (Z.GT . (N#0.005) )G0  TO  1100 

440 

C 

441 

C 

DOUBLE  STEP  SIZE 

442 

C 

443 

DX=2.0*DX 

444 

DO  1070  J2=l ,  J 

445 

GH( J2)=GH( J2)*GH( J2) 

446 

1070 

CONTINUE 

447 

DO  1090  J3=1,J 

448 

DO  1080  J4=1,J 

449 

XI 1(  J3,  J4)=X11( J3,J4)*X11 ( J3,J4) 

450 

X22« J3,J4)=X22( J3,J4)*X22( J3, J4) 

451 

1080 

CONTINUE 

452 

1090 

CONTINUE 

453 

PRINT  1106 , DX 

454 

1106 

FORMAT ( 1X,3HBX=»D17. 10) 

455 

1100 

RN=R 

456 

IF (RMAX-(R-RO) >1110,590,590 

457 

1110 

STOP 

458 

END 

APPENDIX  B 

SAMPLE  OUTPUT  FROM  FAW 


FREQUENCY3  0 . 3000000000D  03 
SOUND  SPEED=  0 . 1500000000D  04 
BETA3  0 . 3500000000D  01 

NORMALIZATION  CONST  ANT  =  0 . 6670000000D-01 

INITIAL  STEP  SIZE3  0 . 1 OOOOOOOCOD  02 

MAXIMUM  DISTANCE=  0 . 1000000000D  05 

SOURCE  SIZE=  O.OOOOOOOOOOD  00 

SPREADING  FACTOR=  0.00 

ABSORPTION  FLAG=  0 

NUMBER  OF  HARMONICS=  40 

INITIAL  G  COEFFIC IENTS=  1 

INITIAL  H  COEFFICIENTS3  0 

NUMBER  OF  HARMONICS  OUTPUT=  5 

PRINT  OUT  INTERVAL3  100 

ALPHA ( 1 ) =  0. 1120000000D-05 

NUMBER  OF  MODIFIED  ABSORPTION  COEFFICIENTS= 
INPUT  WAVEFORM 
G(  1>=  0. 6670000000 D-01 

BIST ANCE3  0. 1000000000D  04 
SIN  COS 

1  0 . 9941231 144D  00  O.OOOOOOOOOOD  00 

2  0.9622622120 D-01  O.OOOOOOOOOOD  00 

3  0. 1394984526D-01  O.OOOOOOOOOOD  00 

4  0.2394449333D-02  O.OOOOOOOOOOD  00 

5  0.451 1965196 D -03  O.OOOOOOOOOOD  00 

DISTANCE3  0. 2000000000 D  04 

SIN  COS 

1  0.9788882494D  00  O.OOOOOOOOOOD  00 

2  0. 1845999987D  00  O.OOOOOOOOOOD  00 

3  0. 519703909411-01  O.OOOOOOOOOOD  00 

4  0. 1729678606D-01  O.OOOOOOOOOOD  00 

5  0.6313700593D-02  O.OOOOOOOOOOD  00 

DISTANCE3  0.30000000000  04 

SIN  COS 

1  0. 9546672934 D  00  O.OOOOOOOOOOD  00 

2  0 . 2587699767D  00  O.OOOOOOOOOOD  00 

3  0. 1041698315D  00  O.OOOOOOOOOOD  00 

4  0. 4945203332 D-Oi  O.OOOOOOOOOOD  00 

5  0.257099491  ID -01  O.OOOOOOOOOOD  00 

DX=  0 . 2000000000D  02 

DISTANCE3  0. 4000000000 D  04 


SIN 

COS 

1 

0 . 92202 1 4992D 

00 

O.OOOOOOOOOOD 

00 

n 

0.3138414940D 

00 

O.OOOOOOOOOOD 

00 

3 

0. 1574155764D 

00 

O.OOOOOOOOOOD 

00 

4 
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